# This file contains the S-plus code that was used to implement the estimators and functionals in "Local linear variable bandwidth hazard rate estimation", Journal of Statistical Computation and Simulation, (81): 1517-1531, (2011)
# As a fully working example the code herein can estimate the hazard rate from the chi square distribution. 
# More general use of the estimates is also feasible by supplying the corresponding data set.  Use of the functions is self - explanatory and is further facilitated by the comments at he preamble of each function.


"LambdaHat"<-
#Watson & Leadbetter estimator (Hazard Analysis I, Biometrika (1964))  -order statistics version. 
function(xin, xout, h, kfun)
{
    n<- length(xin)
    n1<-1:n
    xin.use<-sort(xin)
    arg1<-(sapply(xout, "-", xin.use))/h	
    arg2<-kfun(arg1)* 1/(n-n1+1)
    arg3<-colSums(arg2)* 1/h
    arg3
}

"Epanechnikov"<-   #Epanechikov kernel, square root of 5 is taken as 2.236068 (to save comp. time).
function(x){ ifelse(abs(x)< 2.236068, (3/4)*(1-((x^2)/5 ))/2.236068, ifelse(abs(x)>2.236068, 0,0)) }

"IntEpanechnikov"<-function(x) { ifelse(x< -2.236068, 0, ifelse(x> 2.236068, 1, .5- (2.236068*x*(x^2-15))/100)) }


"IntKde"<-function(xin, xout, h, kfun)
{
	n<- length(xin)
	arg1<-(sapply(xout, "-", xin))/h	
	arg2<-kfun(arg1)
	arg3<-colSums(arg2)
	arg3/n
}
	
"PilotWidth1"<- function(sample)
{
    #returns Silverman's default width 1.06 A n^(-1/5)
    n<-length(sample)
    StD<-stdev(sample)
    IqR<-diff(quantile(sample, c(.25, .75)))
    tmp<-min(StD, IqR/1.34)
    DefWidth<- 1.06*tmp*n^(-1/5)
    DefWidth
}

"sn.0"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun) 
											 sum(1/h *kfun((xin-xout[i])/h)), xin, xout, h, kfun) }

"sn.1"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun)  
											 sum( 1/h * kfun((xin-xout[i])/h)*(xin-xout[i])), xin, xout, h, kfun) }

"sn.2"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun) 
											 sum(1/h *kfun((xin-xout[i])/h)*((xin-xout[i])^2)), xin, xout, h, kfun) }


"tn.0"<-function(xin, xout, h, kfun, Y) { sapply(1:length(xout), function(i, xin, xout, h, kfun,Y)  
                                          sum(1/h *kfun((xin-xout[i])/h)*Y), xin, xout, h, kfun, Y) }

"tn.1"<-function(xin, xout, h, kfun, Y){sapply(1:length(xout), function(i, xin, xout, h, kfun,Y)  
                                        sum(1/h *kfun((xin-xout[i])/h)*((xin-xout[i])*Y)),xin, xout, h, kfun, Y) }


"HazardEstimate"<-function(BinCenters, xout, BandUse, kfun, ci)
{
	sn0<-sn.0(BinCenters, xout, BandUse, kfun)
	sn1<-sn.1(BinCenters, xout, BandUse, kfun)
	sn2<-sn.2(BinCenters, xout, BandUse, kfun)
	tn0<-tn.0(BinCenters, xout, BandUse, kfun, ci)
	tn1<-tn.1(BinCenters, xout, BandUse, kfun, ci)
	estuse<-round((tn1*sn1-tn0*sn2)/(sn1^2-sn0*sn2), digits =4)
	estuse
}

"LocLinVBEst"<-function(xin, xout, h, kfun, IntKfun)
{
	MinX<-min(xin)
	MaxX<-max(xin)
	N<-length(xin)
	n<-length(xout)
	Delta<-(MaxX-MinX)/n
	PartInt<-c(MinX, MinX+(1:n)*Delta) #Partion the interval (x[1], x[n]) into subintervals
	FirstBinCenter<-(PartInt[1]+PartInt[2])/2
	BinCenters<- c(FirstBinCenter, FirstBinCenter+(1:(n-1))*Delta)
	fi<-sapply(1:n, function(i, xin, PartInt) length(which(xin > PartInt[i] & xin <= PartInt[i+1])), xin, PartInt) #crete the bin counts
	citmp<-fi/(N-cumsum(fi)+1)
	ci<-citmp/Delta
	lambdahat <- LambdaHat(xin, xin, .6, kfun)   # h instead of .6
	sqrtlhat <-  lambdahat^{1/2} 
	BandUse<- sqrtlhat #.6	
	llest<-HazardEstimate(BinCenters, xout, BandUse, kfun, ci)
	llest
}

#####################################ORDINARY LOCAL LINEAR ESTIMATE CODE#################################

"Ldoubledash"<-function(x) { dnorm(x,0,1)*(x^2-1) }

"PsiEst"<-function(x,g)
#estimate of the \psi_2 functional to be used for the altman method.
{
	n<-length(x)
	Muse<-sapply(xin, "-", xin)	
	sum(Ldoubledash(Muse/g))/(n^2)
}
"sn.01"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun) sum(kfun((xin-xout[i])/h)), xin, xout, h, kfun) }

"sn.11"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun)  sum(kfun((xin-xout[i])/h)*(xin-xout[i])), xin, xout, h, kfun) }

"sn.21"<-function(xin, xout, h, kfun) { sapply(1:length(xout), function(i, xin, xout, h, kfun) sum(kfun((xin-xout[i])/h)*((xin-xout[i])^2)), xin, xout, h, kfun) }

"sn.31"<-function(xin, xout, h, kfun){sapply(1:length(xout), function(i, xin, xout, h, kfun)  sum(kfun((xin-xout[i])/h)*((xin-xout[i])^3)), xin, xout, h, kfun) }

"sn.41"<-function(xin, xout, h, kfun) {sapply(1:length(xout), function(i, xin, xout, h, kfun)  sum(kfun((xin-xout[i])/h)*((xin-xout[i])^4)), xin, xout, h, kfun) }

"tn.01"<-function(xin, xout, h, kfun, Y) { sapply(1:length(xout), function(i, xin, xout, h, kfun,Y)  sum(kfun((xin-xout[i])/h)*Y), xin, xout, h, kfun, Y) }

"tn.11"<-function(xin, xout, h, kfun, Y){sapply(1:length(xout), function(i, xin, xout, h, kfun,Y)  sum(kfun((xin-xout[i])/h)*((xin-xout[i])*Y)),xin, xout, h, kfun, Y) }

"tn.21"<-function(xin, xout, h, kfun, Y){sapply(1:length(xout), function(i, xin, xout, h, kfun,Y) sum(kfun((xin-xout[i])/h)*((xin - xout[i])^2 )*Y),xin, xout, h, kfun, Y) }


"DerivativeEstimate"<-function(BinCenters, xout, h, kfun, ci)
{
	sn0<-sn.01(BinCenters, xout, h, kfun)
	sn1<-sn.11(BinCenters, xout, h, kfun)
	sn2<-sn.21(BinCenters, xout, h, kfun)
	sn3<-sn.31(BinCenters, xout, h, kfun)
	sn4<-sn.41(BinCenters, xout, h, kfun)
	tn0<-tn.01(BinCenters, xout, h, kfun, ci)
	tn1<-tn.11(BinCenters, xout, h, kfun, ci)
	tn2<-tn.21(BinCenters, xout, h, kfun, ci)
	top<- tn2*((sn1^2)-sn2*sn0)+tn1*(sn3*sn0-sn1*sn2)+tn0*((sn2^2)-sn1*sn3)
	bot<-(sn3^2) * sn0 -2*(sn1 *sn2 * sn3) - (sn0*sn2*sn4) + (sn1^2) *sn4 + sn2^2
	ThetaUse<- 4*sum((top/bot)^2)
	ThetaUse
}
"HazardEstimate1"<-function(BinCenters, xout, BandUse, kfun, ci)
{
	sn0<-sn.01(BinCenters, xout, BandUse, kfun)
	sn1<-sn.11(BinCenters, xout, BandUse, kfun)
	sn2<-sn.21(BinCenters, xout, BandUse, kfun)
	tn0<-tn.01(BinCenters, xout, BandUse, kfun, ci)
	tn1<-tn.11(BinCenters, xout, BandUse, kfun, ci)
	estuse<-round((tn1*sn1-tn0*sn2)/(sn1^2-sn0*sn2), digits =4)
	estuse
}


"AltmanBand"<-function(x,g) { (45/(-7*PsiEst(x,g)*length(x)))^(1/3) }


"LocLinEst"<-function(xin, xout, h, kfun, IntKfun)
{
	MinX<-min(xin)
	MaxX<-max(xin)
	N<-length(xin)
	n<-length(xout)
	Delta<-(MaxX-MinX)/n
	PartInt<-c(MinX, MinX+(1:n)*Delta) #Partion the interval (x[1], x[n]) into subintervals
	FirstBinCenter<-(PartInt[1]+PartInt[2])/2
	BinCenters<- c(FirstBinCenter, FirstBinCenter+(1:(n-1))*Delta)
	fi<-sapply(1:n, function(i, xin, PartInt) length(which(xin > PartInt[i] & xin <= PartInt[i+1])), xin, PartInt) #crete the bin counts
	citmp<-fi/(N-cumsum(fi)+1)
	ci<-citmp/Delta
	#b<-AltmanBand(xin, (length(xin))^(-.3))
	#Muse<-1/(1-IntKde(xin, max(xout), b, IntKfun)) -1
	#ThetaUse<-DerivativeEstimate(BinCenters, xout, h, kfun, ci)
	BandUse<- h#( (.67*Muse) /( N * .25^2* ThetaUse))^(0.2)
	#cat(BandUse)
	llest<-HazardEstimate1(BinCenters, xout, BandUse, kfun, ci)
 llest  
}


#############################END ORDINARY LOCAL LINEAR ESTIMATE CODE##################################
set.seed(154) 

Iterations<-20
matLVB<-matrix(nrow=80, ncol=Iterations)
matLL<-matrix(nrow=80, ncol=Iterations)

for(i in 1:20)
{
xin<-rchisq(400,12)
xout<-seq(min(xin), 20, length=80)

llvdhre<-LocLinVBEst(xin,xout, 2, Epanechnikov, IntEpanechnikov) #.05
ll<-LocLinEst(xin,xout, 2, Epanechnikov, IntEpanechnikov)
matLVB[,i]<-llvdhre
matLL[,i]<-ll
}

plot(xout, dchisq(xout,12)/(1-pchisq(xout,12)), type="l", ylim=c(0, .35), xlab="Hazard rate", ylab="x")
lines(xout, rowMeans(matLVB), lty=3, col=2, type="l")
lines(xout, rowMeans(matLL), lty=4, col=3, type="l")